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Abstract 

Given a compact semi-algebraic set ScK n and a polynomial map / : R n —>► R m , 
we consider the problem of approximating the image set F = /(S) C R m . This in¬ 
cludes in particular the projection of S on R m for n > m. Assuming that F C B, 
with B C R m being a “simple” set (e.g. a box or a ball), we provide two methods to 
compute certified outer approximations of F. Method 1 exploits the fact that F can 
be defined with an existential quantifier, while Method 2 computes approximations 
of the support of image measures. The two methods output a sequence of super¬ 
level sets defined with a single polynomial that yield explicit outer approximations 
of F. Finding the coefficients of this polynomial boils down to computing an op¬ 
timal solution of a convex semidefinite program. We provide guarantees of strong 
convergence to F in L\ norm on B, when the degree of the polynomial approxima¬ 
tion tends to infinity. Several examples of applications are provided, together with 
numerical experiments. 


Keywords semi-algebraic sets; semidefinite programming; moment relaxations; poly¬ 
nomial sums of squares. 


1 Introduction 

Consider a polynomial map / : R n —> M m , x i-> /(x) := (/i(x),... ,/ m (x)) e M m [x] of 
degree d := max{deg/i,... ,deg/ m } and a compact basic semi-algebraic set 

S := {x G R" : fif(x) > 0,... s(x) > 0} (1) 

defined by polynomials gf ,..., g^ s e M[x]. 
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Since S is compact, the image set 


F := /(S) 

is included in a basic compact semi-algebraic set B, assumed to be “simple” (e.g. a box 
or a ball) and described by 

B := {y € R" : gf (y) > 0,... l9 ®„(y) > 0} (2) 

for some polynomials gf ,... , G K[y], 

The purpose of this paper is to approximate F, the image of S under the polynomial 
map /, with super level sets of single polynomials of fixed degrees. One expects the 
approximation to be tractable , i.e. to be able to control the degree of the polynomials 
used to define the approximations. This appears to be quite a challenging problem since 
the polynomial map / and the set S can be both complicated. 

This problem includes two important special cases. The first problem is to approximate 
the projection of S on M' m for n > m. The second problem is the approximation of Pareto 
curves in the context of multicriteria optimization. In [ Mlll.l 1 . we reformulate this 
second problem through parametric polynomial optimization, which can be solved using 
a hierarchy of semidefinite approximations. The present work proposes an alternative 
solution via approximations of polynomial images of semi-algebraic sets. 

In the case of semi-algebraic set projections, notice that computer algebra algorithms pro¬ 
vide an exact description of the projection. These algorithms are based on real quantifier 
elimination (see e.g. |Tar511 !Col74 . BPR.96 ]). For state-of-the-art computer algebra algo¬ 
rithms for quantifier elimination, we refer the interested reader to the survey [ Ba,sl4l | and 
the references therein. Quantifier elimination can be performed with the famous cylindri¬ 
cal algebraic decomposition algorithm. For a finite set of s polynomials in n variables, the 
(time) computational complexity of this algorithm is bounded by ( sd) 2 ° ( "\ thus doubly 
exponential |Col741 |W76| . In [GJ88 ] an algorithm was proposed to find real elements of 
semi-algebraic sets in sub-exponential time. The Block Elimination Algorithm is a singly 
exponential algorithm to eliminate one block of n variables out of n + rn variables, with 
a complexity bounded by s n + 1 d°( n + m ) ( se e [BPR061 Chapter 14] for the formalization of 
this algorithm). For applications that satisfy certain additional assumptions (e.g. rad- 
icality, equidimensionality, etc.), one can use the variant quantifier elimination method 
proposed in [ HD 12 ], which is less computationally demanding. 

Providing approximation algorithms for quantifier elimination is interesting on its own 
because it may provide simpler answers than exact methods, with a more reasonable 
computational cost. On the one hand, we do not require an exact description of the 
projection but rather a hierarchy of outer approximations with a guarantee of convergence 
to the exact projection. On the other hand, the present methodology only requires the 
following assumptions: 1) the set S is compact and 2) either the semi-algebraic set S or 
F (resp. B\F) has nonempty interior. 


Contribution and general methodology We provide two methods to approximate 
the image of semi-algebraic sets under polynomial applications. 
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• Method 1 consists of rewriting F as a set defined with an existential quantifier. Then, 
one can outer approximate F as closely as desired with a hierarchy of superlevel sets 
of the form F* := {y G B : q r (y) > 0} for some polynomials q r G R[y] of increasing 
degrees 2r. 

• Method 2 consists of building a hierarchy of relaxations for the infinite dimensional 
moment problem whose optimal value is the volume of F and whose optimum is the 
restriction of the Lebesgue measure on F. Then, one can outer approximate F as 
closely as desired with a hierarchy of super level sets of the form Fj := {y G B : 
w r (y) > 1}, for some polynomials w r G R[y] of increasing degrees 2r. 

Method 1 and Method 2 share the following essential features: 

1. The sets F' and ¥ 2 r are described with a single polynomial of degree 2r. 

2. Assuming non-emptiness of the interior of S, resp. of F and B\F, one has 
linv^oo vol(Fj\F) = 0, resp. lim^oo vol(Fj?\F) = 0, where vol(-) stands for the 
volume or Lebesgue measure. 

3. Computing the coefficient vectors of the polynomials (c/ r ) reFI , resp. (u-y) r qi, boils 
down to finding optimal solutions of a hierarchy of semidefinite programs. The 
size of these programs is parametrized by the relaxation order r and depends on the 
number of variables n, the number m of components of the polynomial / as well as its 
degree d. For the hierarchy of semidefinite programs associated with Method 1, the 
number of variables at step r is bounded by ( n+ ™ r +2r ), with (n s +n B + l) semidefinite 

constraints of size at most ( n+r " + ' j. Step r of the semidefinite hierarchy associated 
with Method 2 involves at most + 2 variables, (ns + 1) semidefinite 

constraints of size at most and 2(n B + 1) semidefinite constraints of size at 

most ( m + r ). 

4. Data sparsity can be exploited to reduce the overall computational cost. 

Method 1 relies on the previous study |Lasl5| . in which the author obtains tractable 
approximations of sets defined with existential quantifiers. The present article provides an 
extension of the result of |Lasl5i Theorem 3.4], where one does not require anymore that 
some set has zero Lebesgue measure. This is mandatory to prove the volume convergence 
result. 

In [HLS09], the authors consider the problem of approximating the volume of a general 
compact basic semi-algebraic set. The initial problem is then reformulated as an infinite 
dimensional linear programming (LP) problem, whose unknown is the restriction of the 
Lebesgue measure on the set of interest. The main idea behind Method 2 is a similar 
infinite dimensional LP reformulation of the problem, whose unknown is /q, the restriction 
of the Lebesgue measure on F. One ends up in computing a finite number of moments of 
the measure /io supported on S such that the image of //q under / is precisely /q. Note 
however that there is an important novelty compared with |HLS09| , in which the set under 
study is explicitly described as a basic compact semi-algebraic set (i.e. the intersection of 
superlevel sets of known polynomials), whereas such a description is not known for F. 
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Structure of the paper The paper is organized as follows. Section [2] recalls the basic 
background about polynomial sum of squares approximations, moment and localizing ma¬ 
trices. Section [3] presents our approximation method for existential quantifier elimination 
(Method 1). Section [5] is dedicated to the support of image measures (Method 2). In 
Section [5j we analyze the theoretical complexity of both methods and describe how the 
system sparsity can be exploited. Section [6] presents several examples where Method 1 
and Method 2 are successfully applied. 


2 Notation and Definitions 


Let M[x] (resp. R 2r [x]) be the ring of real polynomials (resp. of degree at most 2r) in the 
variable x = (xi,, x n ) G M". for r G N. With S a basic semi-algebraic set as in (fl]), we 
set rf := \(deggf)/2~\,j = 1 ,... ,n s and with B a basic semi-algebraic set as in (j2J), we 

Let E[x] denote the cone of sum of squares (SOS) 
denote the cone of polynomials SOS of degree at most 2r, 


set rf := \(deggf)/2\,j = l,...,n B 
of polynomial, and let E 


that is E r [x] := E[x] fl R 2r [x] 


For the ease of notation, we set gf(x) := 1 and g B ( y) := 1. For each r G N, let 
Q r (S) (resp. Q r (B)) be the r-truncated quadratic module (a convex cone) generated by 
ffo > • • • > 9ns ( reS P- 9o,---,9nJ- 


Qr(S) 

Qr(B) 



: Si 




: Sj G E r _ r e[y], j = 0,.. 



Now, we introduce additional notations which are required for Method 1. Let us first 
describe the product set 

K := S x B = {(x, y) G R n+m : <?i(x, y) > 0,..., g n k(x, y) > 0} C R n+m , (3) 
with n K := n s + n B and the polynomials gj G M[x, y], j = 1,..., n K are defined by: 

A '' (y) if n s + 1 < j < n K . 

As previously, we set := |’(deg 5 r j)/ 2 ], j = 1,... ,n K , g 0 (x.,y) := 1 and Q r (K) stands 
for the r-truncated quadratic module generated by the polynomials g 0 ,..., g n k. 

To guarantee the theoretical convergence of our two methods, we need to assume the 
existence of the following algebraic certificates of boundedness of the sets S and B: 

Assumption 2.1. There exists an integer j s (resp. j B ) such that gj s = g s := N s — ||x||| 
(resp. g B B = g B := N B — ||y|| 2 ) for large enough positive integers N s and N B . 
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For every a G N n the notation x“ stands for the monomial .i:“' ... x“ n and for every r € N, 
let N” := {« £ N" : J2j =1 a j < r}, whose cardinality is One writes a polynomial 

p£R[x,y] as follows: 

(x, y) p(x, y) = P<*P x “y 13 , 

(a,/9)e N n + m 

and we identify p with its vector of coefficients p = (p a p) in the canonical basis (x“y /3 ), 
a G N n , f3 G N m . 

Given a real sequence z = (z a p), we define the multivariate linear functional 4 : M[x.y] — >■ 
E by 4(p) := Y^apPapZap, for all p G R[x,y], 


Moment matrix The moment matrix associated with a sequence z = (z ai 3)( a , / 3)(=N n + m > 
is the real symmetric matrix M r (z) with rows and columns indexed by N” +m , and whose 
entries are defined by: 

M r (z)((a,/5),(h, 7 )):=4(x“ + V +7 ), Va,5e^, V/3, 7 eN™. 


Localizing matrix The localizing matrix associated with a sequence z = (z a p) ( a ,p)eN n + m 
and a polynomial q G R[x,y] (with g(x, y) = Yl u ,v QuvX u y v ) is the real symmetric matrix 
M r (gz) with rows and columns indexed by N” +m , and whose entries are defined by: 

M r ( (7 z)((a,/5),(h, 7 )):=4(g(x,y)x“ + V +7 ), Vcqh G N?, V/3, 7 eN™. 


We define the restriction of the Lebesgue measure on a subset A C B by \A{dy) := 
1a (y) dy, with 1a : B —> {0,1} denoting the indicator function on A: 


My) 


fl if ye A, 

) 0 otherwise. 


The moments of the Lebesgue measure on B are denoted by 

2 J:=|y%(dy)GE, f3 G hT (4) 

We assume that the bounding set B is “simple” in the following sense: 

Assumption 2.2. The moments Q of the Lebesgue measure on B can be explicitly 
computed using cubature formula for integration. 


3 Method 1: existential quantifier elimination 

3.1 Semi-algebraic sets defined with existential quantifiers 

The set F = /( S) is the image of the compact semi-algebraic set S under the polynomial 
map / : S —* B, thus it can be defined with an existential quantifier: 

F = {yGB:3xGS s.t. h/(x, y) > 0}, 
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with 


hf : R n+m —» R, (x,y)^/i / (x,y) 


-||y- /(x)||a 


m 


-E(% -/j( x )) 2 - 

i=i 


Let us also define 

h : M m ->• R, y i-> /i(y) := sup /i/(x, y). 

xGS 

Theorem 3.1. There exists a sequence of polynomials (p r ) rG N C R[y] such that p r ( y) > 
h/(x, y) /or all r G N, x G S, y G B and such that 

Urn J \p r (y) ~ My) I A B (dy) = 0. (5) 

Proof. The result follows readily from [ Las 15 . Theorem 3.1 (3.4)] with the notations 
x i — y, y i — x, K i — B x S, K x <— S 0 and Jf G- —h, which is lower semi- 
continuous. □ 

Theorem 3.2. For each r G N, define F r := {y G B : p r (y) > 0} where the sequence of 
polynomials ( p r )ren C M[y] is as in Theorem \3.1\ Then F r D F and one has 

lim vol(F r \F) = 0. (6) 

r—¥ oo 

Proof. Let r e N. By assumption, one has p r (y) > h/(x,y) for all x e S,y E B. Thus, 
one has p r ( y) > h(y) for all y G B, which implies that F r D F. 

It remains to prove ©• Let us define F (k) := {y G B : h(y) > —1/k}. First, we show 
that 


lim volF(fc) = volF. (7) 

k—> oo 

For each k G N, one has F( k + 1) C F (k) C B, thus the sequence of indicator functions 
(lF(fc))fceN is non-increasing and bounded. Next, let us show that for all y G B, l F (fe)(y) —> 
l F (y), as k —> oo: 

• Let y G F. By the inclusion F C F (k), for each k G N, l F (fc)(y) = l F (y) = 1 and 
the result trivially holds. 

• Let y G B\F, so there exists e > 0 such that h(y) = —e. Thus, there exists k 0 G N 
such that for all k > k 0 , y G B\F(fc). 

Hence, l F (/t)(y) —> l F (y) for each y G B, as k —> oo (monotone non-increasing). By the 
Monotone Convergence Theorem, l F ( fc )(y) —> l F (y) for the L\ norm on B and ([7]) holds. 

Next, we prove that for each k G N, 

lim volF r < volF(fc). (8) 
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By Theorem |3.1 


0. Thus, by [As 
for every e > 0, 


applied to the sequence (p r )re n, one has lim^oo / |p r (y) — ^(y)| ^B(dy) = 
i72 . Theorem 2.5.1], the sequence converges to h in measure, i.e. 


hm vol({y G B : |p r (y) - h(y)| > e}) = 0. 


(9) 


For every k > 1, observe that: 


volF r = vol(F r n {y G B : |p r (y) - h(y)| > 1/fc}) + 
vol(F r n {y G B : |p r (y) - h(y)\ < 1 /k}). 


( 10 ) 


It follows from ([9]) that lim^oo vol(F r fl {y G B : |p r (y) — h(y)| > 1/A;}) = 0. In addition, 
for all r G N, 

vol(F r n {y G B : |p r (y) - h( y)| < 1/fc}) < vol({y G B : h( y) > —1/fc}) = volF(fc). 

( 11 ) 


Using both (10) and and letting r —>• oo, yields ft Thus, we have the following 
inequalities: 

volF < hm volF r < volF(fc). 

r—OO 

Using ([7]) and letting k —>■ 00 yields the desired result. □ 


3.2 Practical computation using semidefinite programming 


In this section we show how the sequence of polynomials of Theorems |3.1| and |3.2| can 
be computed in practice. Define := max{d, ri,, r n k}. For r > rA J n , consider the 
following hierarchy of semidefinite programs: 


Pr : = in f Z 

s.t. q — hf G Q r (K), 
q G R 2r [y]. 


The semidefinite program dual of © reads: 

d* r := sup £ x (hf) 

Z 

s.t. M r (z) >: 0, 


h 0, j = 1,..., 
4(y = V/3gN™. 


K 


n 


Theorem 3.3. Let r > r[^( u and suppose that Assumption 


2.1 


holds. Then: 


( 12 ) 


(13) 


1. p* = d*, i.e. there is no duality gap between the semidefinite program (12) and its 


dual (13). 
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2 . 


The semidefinite program (13) has an optimal solution, 
nonempty interior, then the semidefinite program (12) has 
and the sequence ( q r )rm converges to h in L\ norm on B; 


In addition, if S has 
an optimal solution q r , 


Um J \q r (y) - h( y)| A B (dy) = 0. 


(14) 


3. Defining the set 

K := {y G B : q r { y) > 0} 

it holds that 

fJdf 

and 

lim vol(F*\F) = 0. 

r»—Vmr'i \ I \ / 


Proof. 1. Let V r (resp. V*) stand for the feasible (resp. optimal) solution set of the 
semidefinite program (jl3j). First, we prove that V r 0. Let z = (z Qj 8)( a)j9 ) eN ™+m 
be the sequence moments of Ak, the Lebesgue measure on K = S x B . Since the 
measure is supported on K, the semidefinite constraints M r (z) y 0, M r _ r K(<yyz) A 

0, j = 1,..., n K are satisfied. By construction, the marginal of Ak on B is A B and 
the following equality constraints are satisfied: £ z (y^) = zf for all ft G Thus, 
the finite sequence z lies in T> r 0. 


Note that Assumption 2.1 implies that the semidefinite constraints M r _i(g s z) y 0 
and M r _ 1 (ry B z) y 0 both hold. Thus, the first diagonal elements of M r _ 1 (r/ S z) and 
M r _i(g B z) are nonnegative, and since £ z (l) = z B , it follows that £ z (x? k ) < ( N s ) k z B , 
i = 1,..., n and I z (yj k ) < ( N B ) k z B , j = 1,..., m, k = 0,..., r and we deduce 
from |LN07l Lemma 4.3, p. Ill] that \z a g\ is bounded for all (a, (3) G Thus, 

the feasible set V r is compact as closed and bounded. Hence, the set 'D* is nonempty 
and bounded. The claim then follows from the sufficient condition of strong duality 
in [TmQijj. 


Assume that S has nonempty interior, so K = S x B has also nonempty interior. 
Thus, the feasible solution z (defined above) satisfies M r (z) y 0, M r _ r .K((yfjz) >- 0, 

j = 1,... ,n K , which implies that Slater’s condition holds for (13). Note also that 
the semidefinite program (12) has the trivial feasible solution <7 = 0 since —hf is 
SOS by construction. As a consequence of a now standard result of duality in 


semidefinite programming (see e.g. [ VB94 ]). the semidefinite program (12) has an 
optimal solution q r e R 2r [y]. 

Let us consider a sequence of polynomials (pk)ke n C K[y] as in Theorem 3.1 Now, 
fix e > 0. By Theorem |3.1[ there exists k 0 G N such that 


J \Pk(y) ~ h(y)\X B (dy) < e/2, 


(15) 


for all k > k 0 . Then, observe that the polynomial p\ := Pk + e/(2volB) satisfies 
p%(y) — hf(x , y) > 0, for all x G S, y G B. For r G N large enough, as a consequence 
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of Putinar’s Positivstellensatz (e.g. [ Las09i Section 2.5]), there exist so,• • •, s„k £ 
E[x, y] such that 

n K 

pK y) - M x > y) = 5>i( x >y)^( x ,y), 

i=o 


with deg(sjgfj) < 2r for j = 0,..., n K . And so, p' k — hf lies in the r-truncated 
quadratic module Q r (K), which implies that p' k is a feasible solution for Prob¬ 
lem (12). Hence, q r being an optimal solution of problem (12), the following holds: 


J 4 r(y)A B (dy) < J Pfc(y)A B (dy) = J [p fc (y) + e/(2vol(B))]A(dy). (16) 


Combining (15) and (16) yields J |g r (y) — h(y)|A(dy) < e, concludes the proof. 


3. This is a consequence of Theorem 3.2 


□ 


4 Method 2: support of image measures 

Given a compact set Ad", let Ad (A) stand for the vector space of finite signed Borel 
measures supported on A, understood as functions from the Borel sigma algebra 13(A) to 
the real numbers. Let C(A) stand for the space of continuous functions on A, equipped 
with the sup-norm (a Banach space). Since A is compact, the topological dual (i.e. the set 
of continuous linear functionals) of C(A) (equipped with the sup-norm), denoted by C( A)', 
is (isometrically isomorphically identified with) Ad (A) equipped with the total variation 
norm, denoted by || • ||tv- The cone of non-negative elements of C(A), resp. Ad (A), 
is denoted by C+(A), resp. Ad+(A). The topology in C+(A) is the strong topology of 
uniform convergence while the topology in Ad+(A) is the weak-star topology (see [Bar02. 
Chapter IV] or [ Lue97 . Section 5.10] for more background on weak-star topology). 

Recall that A B stand for the Lebesgue measure on B. If p, is e Ad (A), the notation 
p, <C v stands for // being absolutely continuous w.r.t. u, whereas the notation ji < u 
means that u — p G Ad+(A). For background on functional analysis and measure spaces 
see e.g. [ RF10 . Section 21.5]. 

4.1 LP primal-dual conic formulation 

Given a polynomial application / : S —> B, the pushforward or image map 

f# : Ad(S) —>• Ad(B) 


is defined such that 

f#p o(A) := /i 0 ({x e S : /(x) e A}) 


for every set A e B(B) and every measure po E Ad(S). The measure fpPo E Ad(B) is 
then called the image measure of p 0 under /, see e.g. |AFP0f)l Section 1.5]. 
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To approximate the image set F = /(S), one considers the infinite-dimensional linear 
programming (LP) problem: 

p* := sup / pi 

w,/u,Ai 9 

s.t. Pi + Pi — A b , (17) 

hi 0) 

/r 0 eAd + (S), e X+(B), 

In the above LP, by definition of the image measure, /ii exists whenever p 0 is given. The 
following result gives conditions for p 0 to exist whenever pi is given. 

Lemma 4.1. Given a measure p, i G Ad + (B), there is a measure p 0 G _A4+(S) such that 
f#Po = Pi if and only if there is no continuous function v G C(B) such that n(/(x)) > 0 
for all x G S and f v(y)dpi(y) < 0. 

Proof. This follows from [ CK77 . Theorem 6] which is an extension to locally convex topo¬ 
logical spaces of the celebrated Farkas Lemma in finite-dimensional linear optimization. 
One has just to verify that the image cone f#(A4 + ( S)) = {f#po ■ Po G A1+(S)} is closed 
in the weak-star topology <r(Al(B),C(B)) of At (B). This. in turn, follows from continuity 
of / and compactness of S. □ 

Lemma 4.2. LP 0 admits an optimal solution (p$, p\, pi). Moreover, p\ = Ap and 
the optimal value of LP ( [17] ) is p* = volF. 

Proof. First, we prove that for p\ = Ap G A4+(B), there is a measure pf G A4 + (S) such 
that /#ho = hi- Indeed, by Radon-Nikodym there exists a function q\ G Li(Ap) such 

that dp\( y) = <?i(y)dAp(y) and q\ (y) > 0 for all y G F. The claim follows then from 

Lemma [4.1| since it is impossible to find a function v G C(B) such that v(y) > 0 for all 
y G F while satisfying f n(y)gi(y)dAp(y) < 0. Define p\ — Ab — p\, then (p^, p\, p\) 
is admissible for LP Exactly as in the proof of [ HLSfifi . Theorem 3.1], one finally 

shows that this triplet is optimum as well as uniqueness of //], yielding the optimal value 
p* = volF. □ 

Next, we express problem ( [T7| ) as an infinite-dimensional conic problem on appropriate 
vector spaces. By construction, a feasible solution of problem 0 satisfies: 

/ v(y)pi(dy) - f u(/(x))ho(dx) = 0, (18) 

J B J S 

J B w(y) Pi(dy) + w( y) pi(dy) = w( y) A(dy), (19) 

for all continuous test functions v,w G C(B). 

Then, we cast problem 0 as a particular instance of a primal LP in the canonical form 
given in [Bar02l 7.1.1]: 

p* = sup ( x,c)i 

X 

s.t. Ax = b, (20) 

x G Ef, 

with 
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• the vector space E\ := A'f(S) x .M(B) 2 ; 

• the vector space F\ := C(S) x C(B) 2 ; 

• the duality (•, -)i : E\ x F\ —> R, given by the integration of continuous functions 
against Borel measures, since E\ = F[\ 

• the decision variable x := (/i 0 , /ii, Ai) £ -^i and the reward c := (0,1, 0) G Up 

• U 2 : = .M(B) 2 , F 2 := C(B ) 2 and the right hand side vector b := (0, A) G E 2 — F' 2 \ 

• the linear operator A : E\ —> E 2 given by 


Al) 


~f#F o + /^i 

Ah + Ai 


Notice that all spaces Ei, E 2 (resp. Pj, P 2 ) are equipped with the weak topologies 
a(Ei,Fi), cf(E 2 , F 2 ) (resp. a(Fi,Ei), a(F 2 ,E 2 )). Importantly, a(Ei,Fi) is the weak-star 
topology (since E\ = F[). Observe that A is continuous with respect to the weak topology, 
as A'(F 2 ) C F i. 

With these notations, the dual LP in the canonical form given in [ Bar02 . 7.1.2] reads: 


d* = inf ( b , y) 2 

s.t. A' y — c e C+(B) 2 


( 21 ) 


with 


• the dual variable y := (v,w) G E 2 ; 

• the (pre)-dual cone C + (B) 2 , whose dual is E j 1 "; 

• the duality pairing (•, -) 2 : E 2 x F 2 —>■ R, with E 2 = F 2 , 

• the adjoint linear operator A! : F 2 —y F\ given by 


A! (v , w) 


-vo f 
V + w 
w 


Using our original notations, the dual LP of problem ( fl7| ) then reads: 

d* : =inj j w(y) X B (dy) 

s.t. n(/(x)) >0, Vx G S, 

w(y) > 1 + v(y), Vy G B, 
w(y) >0, Vy G B, 
v,w G C(B). 


( 22 ) 


11 








Theorem 4.3. There is no duality gap between problem ( [17] ) and problem (22), i.e. p 
d*. 


Proof. This theorem follows from the “zero duality gap” result from [Bar02 . Theorem 
7.2], if one can prove that the cone 

A(Ef):={(Ax,(x,c) 1 ):xeEf} (23) 

is closed in E 2 x M. To do so, let us consider a sequence (x^) = (pq , p^)) C Ef 

such that Ax^ k > —> s = (si, s 2 ) and (x^ k \c) i —> t. Let us prove that (s, t) = (Ax*, (x*, c)i) 
for some x* = (p^, p\, pf) G Ef. As c = (0,1,0), one has ||/4 ^IItv = Ib Pi (^y) t(> 

0), thus sup fc \\p^ ||tv < oo. Therefore there is a subsequence (denoted by the same 
indices) (p^) which converges to //] G _A4+(B) for the weak-star topology In particular 
\\pI\\tv = t. Hence from pi + pi —>■ s 2 one deduces that p^ — t s 2 — p\ for the weak- 
star topology. But then we also have t Si — p\ in the weak-star topology of 

At(B). Therefore —si + p\ is a positive measure. So let p$ be such that f#Py = —s\ + p\ 
guaranteed to exist since we have seen that f#(Ai+(S)) is weak-star closed. Then we 
have A(x*) = (si,s 2 ) and (x*,c) = t, the desired result. □ 


4.2 Practical computation using semidefinite programming 


For each r > r^ n := max{ \rf/d ],..., \r%/d],rf ,... ,r® B }, let z 0 = (z 0 /?)/3eN™ be the 
finite sequence of moments up to degree 2r of measure p 0 . Similarly, let zi and zi stand 
for the sequences of moments up to degree 2r, respectively associated with pi and p±. 
Problem can be relaxed with the following semidefinite program: 


P* ■= SUp Z 10 

Z 0 ,Zl,Zi 

s.t. z ip + z ip = zf, L Z0 (f(xf) = zip, V/3 G N£, 
M rd-rf(9j z o) SO, j = 0,..., n s , 

Mr-rf (hfzi) y 0, M r _ r B(^Zi) so, j = 0,... ,n B . 


(24) 


Consider also the following semidefinite program, which is a strengthening of problem (22) 
and also the dual of problem (25): 


d * ■= inf E w P z f 

S.t. V O f £ Qrd(S), 

w — 1 — v G Q r (B), 
w G Q r (B), 
v,w e M 2r [y]. 


(25) 


Theorem 4.4. Let r > r^L and suppose that both F and B\F have nonempty interior 
and that Assumption \2.1\ holds. Then: 
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1. p* — d*, i.e. there is no duality gap between the semidefinite program ( 24 ) and its 
dual (25). 


2. The semidefinite program (25) has an optimal solution ( v r ,w r ) G R 2 r[y] x M 2 r[y]? 
and the sequence (w r ) converges to lp in Li norm on B: 


™ / My) _ lr ( y ) AB < <iy ) = 


(26) 


3. Defining the set 
its holds that 

and 


F 2 r := {y G B : w r (y) > 1} 

fJdf 

lim vol(F 2 \F) = 0. 


Proof. 


1. Let Hi = Ap, let no be such that f#Ho — hi as l n Lemma 4.1, and let fii = Ab — hi 


ToT 

so that (/i 0 ,hi,hi) I s feasible for LP (17). Given r > rU' n , let z 0 , Zi and zi be 


the sequences of moments up to degree 2 r of /v, 0 , p,i an( f hi> respectively. Clearly, 
(z 0 ,zi,zi) is feasible for program (24). Then, as in the proof of the first item of 
Theorem |3.3[ the optimal solution set of the program (24) is nonempty and bounded, 
which by |Trnf)5j implies that there is no duality gap between the semidefinite 


program (25) and its dual (24). 


2. Now, one shows that (z 0 ,zi,zi) is strictly feasible for program (24). Using the fact 
that 


(a) F (resp. B\F) has nonempty interior, 

(b) Tii (resp. zi) is the moment sequence of pn (resp. /q), 

one has M r (p^Zi) y 0 (resp. M r (p®zi) y 0), for each j = 0,... ,n B . Moreover, 
M r {gj z 0 ) y 0, for all j = 0,..., n s . Otherwise, assume that there exists a nontrivial 
vector q such that M r (g®zo)q = 0 for some j. As F has nonempty interior, it 
contains an open set A C M m . By continuity of /, the set / -1 (A) := {x G S : 
/(x) G A} is an open set of S and /i 0 (/ -1 (A)) = /q(A) > 0. Then, 


0 = (q, M r (g®z 0 ) q) = f g(x) 2 p s (x) d/r 0 (x) > f q(x) 2 p s (x) d/i 0 (x), 

which yields g(x) 2 gf(x) = 0 on the open set / -1 (A), leading to a contradiction. 

Therefore, as for the proof of the second item of Theorem |3.3[ we conclude that the 
semidefinite program (25) has an optimal solution (v r ,w r ) G R 2 r[y] x ® 2 r[y]- 
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Next, one proves that there exists a sequence of polynomials {w k ) ke n C R[y] such 
that w k (y) > l F (y), for all y G B and such that 

hm J\w k (y) - l F (y)| A B (dy) = 0. (27) 


The set F being closed, the indicator function lp is upper semi-continuous and 
bounded, so there exists a non-increasing sequence of bounded continuous func¬ 
tions h k : B —» R. such that h k ( y) i l F (y), for all y G B, as A: —> oo. Using 
the Monotone Convergence Theorem, h k —> 1 F for the L\ norm. By the Stone- 
Weierstrass Theorem, there exists a sequence of polynomials (w' k ) ke ^ C R[y], such 
that sup yeB \w k (y) — h k ( y)| < 1/k. The polynomial w k := w k + 1/k satisfies 
w k > h k > lp and ( |27| ) holds. 


Let us define w k ■— w k + e/(2 volB), v k := w k — 1. Next, for r e N large enough, one 
proves that (%, is a feasible solution of (25). Using the fact that w k > w k > lp, 
one has w k G Q r (B), as a consequence of Putinar’s Positivstellensatz. For each 
x G S, h fc (/(x)) = w k ( y) - 1 > l F (/(x)) - 1 + e/(2volB) > 0, so v k o / lies 
in Q rt i(S). Similarly, w k — v k — 1 G Q r (B). Then, one concludes using the same 
arguments as for (14) in the proof of the second item of Theorem |3.3| 


3. Let y G F. There exists x G S such that y = /(x). Let ( v r ,w r ) G R 2 r[y] x R 2 r[y] 
be an optimal solution of (25). By feasibility, w r (y) — 1 > v r (y) = iy(/(x)) > 0. 
Thus, F^ D F. Finally, the proof of the convergence in volume is analogous to the 
proof of (|6j) in Theorem 3.2, 


□ 


5 Computational considerations 


5.1 Complexity analysis and lifting strategy 

5.1.1 Method 1 


First, consider the semi definite program (|13|) of Method 1. For r > r^ n , the number of 


variables (resp. size of semidefinite matrices m^) of problem (13) satisfies: 

'n + m + 2 


n 


(!)< 


2 r 


Problem (13) involves (n B + n s + 1) semidefinite constraints of size rn (1> bounded as 
follows: 

+ m + r\ 


m 


(!)< 
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5.1.2 Method 2 


Now, consider the semidefinite program (24). For r > r 1TI / n , the number of variables n (2 ^ 
of Problem (24) satisfies: 


n ^ < 


fn + 2 rd \ fm + 2 r\ 
y 2 rd ) y 2r ) 


Problem ( |24| ) also involves (n s + 1) semidefinite constraints of size at most and 

2(n B + 1) semidefinite constraints of size at most . 

Due to the dependence on the degree d of the polynomial application, one observes that 
the number of variables (resp. constraints) can quickly become large if d is not small. An 
alternative formulation to limit the blowup of these relaxations is obtained by considering 
yi,... , ym as “lifting” variables, respectively associated with /i,..., f m , together with the 
following 2m additional constraints: 

c/n s +j(x, y) := y 3 - /i(x), c/ns+ 2 j( x , y) : = /j( x ) - Vj, j = 


By considering the basic compact semi-algebraic set S]jf t C M n+m given by 
Slift := {( x ,y) e K : ^®S +I (x,y) > 0,... ,^s +2m (x,y) > 0}, 


problem (25) becomes: 


inf 

V,W 

E w p z p 



S.t. 

V G Qr(Sljft), 


w — 1 — V G Q 


W G Q r (B), 


V,w G M 2r [y], 


which is actually equivalent to the following problem: 


inf 


s.t. 


Z w P z f 

pe N™ 

W 1 G Qr(^lift)) 

w e Q r (B), 
w G R 2r [y]. 


(28) 


(29) 


(30) 


The minimal relaxation order of problem (|30| is := max{ f|], rf,..., rf iS , rf,..., r B B } 
and the number of variables n| if ( is bounded as follows: 

n ^ 2 ) < 

n iift ^ 



Problem (30) involves (n s + 2m + 1) semidefinite constraints of size at most y n+ ™ +I J and 

(n B + 1) semidefinite constraints of size at most ( m + r 'j ■ When rn is small and d is large, 
this seems to be a suitable choice to reduce the computational cost of the semidefinite 
program (|24]). Experimental results described further (see Table [l] in Section 6.1 and 
Table [2] in Section 6.4) agree with this observation. 
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5.2 Exploiting sparsity 

As explained above, both Method 1 and Method 2 are computationally demanding in gen¬ 
eral. However, if the polynomials x i->- f 3 (x), (j = 1,,m) have some structured sparsity, 
then one can still exploit sparsity in a way similar to the one described in |WKKM06 . 
ILasOfij to handle problems in higher dimensions. In particular, let {1,..., n} be the union 
Uj=i Ij of subsets Ij C {1,... ,n} and assume that for each j = 1,...,m, the polynomial 
fj involves only variables {ay | i e Ij}. One also suppose that the collection {fi ,..., I m } 
satisfies the so-called running intersection property. 

Definition 5.1. The collection {A,..., I m } of subsets Ij C {1 ,... ,n} satisfies the run¬ 
ning intersection property if the following holds for each j = 1,..., m — 1: 

j 

Ij + 1 0 [J 4 C /j for some l < j. 

k =l 

The following assumption allows to apply the sparse representation result of jhasfifil Corol¬ 
lary 3.9] while using either Method 1 or Method 2. 

Assumption 5.2. The index set {1,... ,n} is partitioned into m disjoint sets Ij, j = 
1,..., m so that: 

1. The collection {I \,..., I m } satisfies the running intersection property. 

2. For each j = 1,... ,n s , there exists some k 3 such that the polynomial gj in |l| 
involves only variables {ay | i 6 I kj }. 

3. In the definition ([Tj) of S, we replace the inequality constraint N s — ||x||| > 0 by 
the m quadratic constraints: 

Nj - X 1 > j = l,..., m. 

ieij 


For each j = 1,..., m, index the variable ijj by n+j and define Ij := Ij U{n+ 1,..., n + 
m}. 


Proposition 5.3. Under Assumption 5.2, the collection {/{' ,... ,1$} of subsets Ij ^ C 
{1,..., n, n + 1,..., n + m} satisfies the running intersection property. 


Proof. The collection {fi ,..., I m } of subsets Ij C {1,..., n} satisfies the running intersec¬ 
tion property. For each j = 1,..., m — 1, there exists l < j such that I :)+ \ fl Ui=i A A h- 

Thus, nui =1 4 1} = (Ij +iCULi h) u{n+l, • • •, n+m} C h U{n+1,..., n+m} = j/ 1} , 
the desired result. □ 


Then Assumption 5.2 allows one to apply the sparse representation result of [ LasOfil 
Corollary 3.9] to the semidefinite program (12) associated with Method 1. Indeed, observe 
that the polynomial (x,y) i-y h^(x,y) can be decomposed as hf = Yffj=\hf-, where for 
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each j = 1,... ,m, the polynomial hf . involves only the variables {ay | i G I j } (the same 
variables involved in fj ) and y. 


Under Assumption |5.2[ this sparse representation result can also be applied for the 
urogram (30) associated with the lifting variant of Method 2 described in 


semidefinite 
Section 


and g s 


5.1.2 


This is due to the fact that for each j = 1 ,... ,m, the polynomials gf tS + . 


7l S +2 j 


involve only the variables {ay | i G I 3 } and y. 


6 Application examples 


Here we present some application examples together with numerical results. In particular, 
this section illustrates that our methodology is a unified framework which can tackle 
important special cases: semi-algebraic set projections (Section |6.2[ ) and Pareto curves 
approximations (Section |6.3 ). Moreover, the framework can be extended to approximate 
images of semi-algebraic sets under semi-algebraic applications (Section |6.4[). 


The numerical results are given after solving either the semidefinite program (12) for 


Method 1, the semidefinite program (25) for Method 2 or the semidefinite program (30) 


for the lifting variant of Method 2, with the Yalmip toolbox |L504] for Matlab. As 


explained in Section 5.1.2 , the outer approximations obtained by Method 2 and its lifting 
variant are the same, but their semidefinite formulations differ. 

Benchmarks are performed on an Intel Core i5 CPU (2.40 GHz) with Yalmip interfaced 
with the semidefinite programming solver Mosek jAAOOj . 


6.1 Polynomial image of semi-algebraic sets 

Example 1. Consider the image of the two-dimensional unit ball S := {x G M 2 : ||x||| < 
1} under the polynomial application /(x) := (x\ + X\X 2 ,X 2 — x\)/2. We choose B = S 
since it can be checked that F = /( S) C B. 


On Figure [I] resp. [2j we represent in light gray the outer approximations F) resp. F 2 
of F obtained by Method 1 resp. 2, for increasing values of the relaxation order r. On 
each figure, the black dots correspond to the image set of the points obtained by uniform 
sampling of S under /. We observe that the approximations behave well around the 
locally convex parts of the boundary of F, and that it is not straightforward to decide 
whether Method 1 or Method 2 provides the best approximations. 


We indicate in Table [l] the data related to the semidefinite programs solved by Mosek 
to compute approximations of increasing degrees, while using Method 1, Method 2 and 
Method 2 with the lifting strategy (see Section 5.1.2 for more details). For each problem, 
“vars” stands for the total number of variables and “size” stands for the size of the 
semidefinite matrices. The computational timings of Method 2 with the lifting strategy 
are similar to those of Method 1, for r = 1,..., 5. However, for r = 6 the size of the 
problem is significantly smaller with Method 1. 
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h h »t h 


(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 

Figure 1: Outer approximations F'. (light gray) of F (black dot samples) for Example [l], 
for r = 1, 2, 3, 4. 



h h h h 


(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 

Figure 2: Outer approximations F^ (light gray) of F (black dot samples) for Example [lj 
for r = 1, 2, 3, 4. 


Table 1: Comparison of timing results for Example |lj 


relaxation order r 

i 

2 

3 

4 

5 

6 


vars 

40 

212 

1039 

4211 

14028 

40251 

Method 1 

size 

30 

111 

350 

915 

1991 

3822 


time (s) 

0.64 

0.72 

0.77 

1.69 

8.22 

40.37 


vars 

286 

2140 

8241 

22720 

51166 

100626 

Method 2 

size 

129 

471 

1029 

1803 

2793 

3999 


time (s) 

0.65 

0.74 

1.54 

3.4 

12.89 

43.74 

Method 2 
with lifting 

vars 

51 

308 

1499 

5882 

19546 

56710 

size 

32 

157 

536 

1411 

3128 

6127 

time (s) 

0.58 

0.66 

0.68 

1.93 

10.07 

63.88 


6.2 Projections of semi-algebraic sets 


For n > m, we focus on the special case of projections. Let / be the projection of S with 
respect to the m first coordinates, i.e. /(x) := (xi,... ,x m ). It turns out that in this 


case, the semidehnite program (30) associated to the lifting variant of Method 2, has the 
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following simpler formulation: 


pe N™ 

w-leQ r (S), (31) 

w e Q r (B), 

W £ E 2 r[^l, • • - 

Example 2. Consider the projection F on the first two coordinates of the semi-algebraic 
set S := {x £ M 3 : ||a:||| < 1, 1/4 — [x\ + 1/2) 2 — < 0, 1/9 — ( x\ — 1/2) 4 — x\ < 0}, 

which belongs to B := {x £ R 2 : \\x\\l < 1}. 


inf 

W 

S.t. 



(a) r = 2 


(b) r = 3 


(c) r = 4 


Figure 3: Outer approximations Fj (light gray) of F (black dot samples) for Example [ 2 J 
for r = 2, 3, 4. 



(a) r = 2 


(b) r = 3 


(c) r = 4 


Figure 4: Outer approximations F 2 (light gray) of F (black dot samples) for Example [ 2 J 
for r = 2,3,4. 


Figure [3] resp. [4] displays approximation of the projection of S on the first two coordinates 
with Method 1 resp. 2. 
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6.3 Approximating Pareto curves 


In MlII. 1 1,. we propose a method to approximate Pareto curves associated with bicri¬ 
teria polynomial optimization problems min xe s{(/i(x), / 2 (x))}. The image space M 2 is 
partially ordered with the positive orthant E+, that is, for every yi,y 2 G E 2 , yi > y 2 
stands for y 2 — yi G E+. A point x is called a weakly Edgeworth-Pareto optimal point , 
when there is no x G S such that f 3 (x) < / 7 (x), j = 1,2. The Pareto curve is the set of 
weakly Edgeworth-Pareto optimal points. For more details on multicriteria optimization, 
we refer the interested reader to [ Jahlfij and the references therein. 

The methodology of [MHL14] consists of reformulating the initial bicriteria optimization 
problem to use a hierarchy of semidefinite approximations for parametric polynomial 
optimization problems. Then, one can apply the framework developed in [ Las 10 i and build 
a hierarchy of semidefinite programs, allowing to approximate as closely as desired the 
Pareto curve. Here we propose to study outer approximations of the set F = (/i(S), / 2 (S)) 
since points along the boundary of a tight outer approximation are expected to be close 
to the Pareto curve. 

Example 3. Let consider the two-dimensional nonlinear problem proposed in [ WCSF01 ]: 
min xe s{(./i(x), / 2 (x))}, with /i(x) := + ( X2 - Xl + 3) 2 , / 2 (x) := x x + x\ and 

S := {x G E 2 : ~( Xl - 2) 3 /2 - x 2 + 2.5 > 0, -x x -x 2 + 8(x 2 - Xl + 0.65) 2 + 3.85 > 0}. 
Instead of /i, we consider /i := (/i(x) — a\)/(hi — ai), where a\ and b\ are given by 
ai := min xe s/i(x) and b\ := /i(x) with x a solution of min xe s/ 2 (x). Similarly, we 
consider a scaled criterion f 2 defined from f 2 . A preprocessing step consists in computing 
lower and upper bounds of the polynomial f\ (resp. f 2 ) over S to define / = (/i,/ 2 ). 
Doing so, one ensures that /(S) is a subset of the unit ball B and the present methodology 
applies. 




>> 




h 

(b) r = 2 


(c) r = 4 


Figure 5: Outer approximations F'. (light gray) of F (black dot samples) for Example |3j 
for r = 1, 2, 4. 

In this case, the Pareto curve is nonconvex and disconnected. As depicted in Figure [5] 
and Figure [6| it is difficult to obtain precise approximations of the whole image set, in 
particular for the subset F Pi B 0 , with an ellipse B 0 := {x G E 2 : (( Xl — 13.7)/1.7) 2 + 
((x 2 — 1.8)/0.5) 2 < 1}. Figure [7] displays more precise outer approximations of degree 8 
(a) and degree 10 (b). 
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(b) r = 2 


(c) r = 4 


Figure 6: Outer approximations (light gray) of F (black dot samples) for Example |3j 
for r = 1, 2, 4. 



Figure 7: Outer approximations Fj (light gray) of F fl B 0 (black dot samples) for Exam¬ 
ple [3j for r = 4, 5. 


6.4 semi-algebraic image of semi-algebraic sets 

Given a semi-algebraic set S as in 0, Methods 1 and 2 can be extended to approx¬ 
imate the image of S under a semi-algebraic application / = (/i,..., f m ). To do so, 
we follow [FPTnj and introduce lifting variables to represent non-polynomial components 
involved ia fi,..., f m , as well as additional polynomial constraints. 

Proceeding as in [LP10 ;. for each semi-algebraic function fj, one introduces additional 
variables x- 7 := (x{,..., xj.) such that the graph {(x, /)(x)) : x e S} = {(pc,x J tj ) : 
(x.x J ) e Sj} for some semi-algebraic set Sj C M n+ b. In the end, one works with the 
lifted set S := {(x, x 1 ,..., x m ) : (x, x J ) e Sj , j — 1,..., m }. 

Example 4. Here, we consider the image of the two-dimensional unit bah S := {x G i 2 : 

11x11 2 < 1} under the semi-algebraic application /(x) := (min(xi + x\X 2 ,x\),X 2 — x\)/2>. 
Remind that 2min(o, b) = a + b — \a — b\, so that 2min(xi + x\X 2 ,x \) = x\ + X\X 2 + 
x\ — | x i + x i X‘i — x\\. To handle the absolute value, we introduce an additional variable 
x 3 together with the equality constraint x\ = (x\ + X\X 2 2 ) 2 and the inequality constraint 
x 3 > 0. 
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-1 - 0.5 0 0.5 1 -1 - 0.5 0 0.5 1 -1 - 0.5 0 0.5 1 -1 - 0.5 0 0.5 1 


h h h h 

(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 

Figure 8: Outer approximations F'. (light gray) of F (black dot samples) for Example [fj 
for r = 1,2, 3, 4. 



h h h h 


(a) r = 1 (b) r = 2 (c) r = 3 (d) r = 4 

Figure 9: Outer approximations F^ (light gray) of F (black dot samples) for Example [fj 
for r = 1,2, 3, 4. 


As for Example [lj we report in Table [2] the data related to the semidehnite problems 
solved by Mosek to compute approximations of increasing degrees, while using Method 
1, Method 2 and Method 2 with the lifting strategy. Method 2 fails to compute polynomial 
approximations of degree higher than eight (r = 4), the system running out of memory 


(indicated with the symbol ”). The lifting strategy described in Section 5.1.2 
this practical limitation. 


overcomes 


Table 2: Comparison of timing results for Example [ 4 ] 


relaxation order r 

1 

2 

3 

4 

5 


vars 

66 

438 

3137 

16993 

73213 

Method 1 

size 

45 

226 

1008 

3387 

9075 


time (s) 

0.68 

0.85 

1.16 

14.41 

147.32 


vars 

715 

12243 

89695 

— 

— 

Method 2 

size 

295 

1957 

6283 

— 

— 


time (s) 

0.83 

3.29 

52.55 

— 

— 


vars 

78 

540 

3788 

20216 

87475 

iviei;noa z 

size 

51 

273 

1262 

4247 

11508 

with lifting 

time (s) 

0.68 

0.96 

1.83 

14.44 

174.80 
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7 Discussion and conclusion 


In this work, we propose two methods to approximate polynomial images of basic com¬ 
pact semi-algebraic sets, a numerical approximation alternative to exact computer algebra 
methods when the latter are too computationally demanding. In its present form, this 
methodology is applicable to problems of modest size, except if some sparsity can be 
taken into account, as explained earlier. Therefore, to handle larger size problems, the 
methodology needs to be adapted. A topic of further investigation is to search for alter¬ 
native positivity certificates, less demanding than the SOS certificates used in this paper 
but more efficient than the LP based certificates as defined in [ Han88 . Vasj . On the one 
hand, the latter are appealing since they yield a hierarchy of LP relaxations (as opposed 
to semidefinite relaxations as in this paper). Moreover, today’s LP solvers can handle 
huge size LP problems, which is far from being the case for semidefinite solvers. On the 
other hand, it has been shown in |Las09j that generically finite convergence cannot occur 
for convex problems, except for the linear case. Finally, it could be interesting to look at 
various compactification procedures to study how the methodology could be generalized 
to non-compact situations. 
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